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SUMMARY 

A finite-element  formulation  of  the  equations  governing  the  laminated 
anisotropic  plate  theory  of  Yang,  Norris  and  Stavsky  is  presented.  The 
theory  is  a generalization  of  Mindlin's  theory  for  isotropic  plates  to 
laminated  anisotropic  plates  and  includes  shear  deformation  and  rotary  in- 
ertia effects.  Finite-element  solutions  are  presented  for  rectangular  plates 
of  antisymmetric  angle-ply  laminates  having  material  properties  that  are 
typical  of  a highly  anisotropic  composite  material.  Two  sets  of  material 
properties  that  are  typical  of  advanced  fiber- reinforced  composites  are  used 
to  show  the  parametric  effects  of  plate  aspect  ratio,  length-to-thickness 
ratio,  number  of  layers,  and  lamination  angle.  The  element  is  also  employed 
to  study  the  bending  of  laminated,  anisotropic  bimodulus-material  plates. 
Results  are  presented  for  single-layer  and  two-layer  cross-ply  rectangular 
plates  subjected  to  sinusoidal  loading. 

The  report  also  presents  a C"  finite  element  for  the  von  Karman  equa- 
tions of  thin  elastic  plates.  The  slope-displacement  relations  are  treated 
as  constraints  using  the  so-called  penalty  method  of  Courant.  The  resulting 
element  contains  the  transverse  deflection  and  two  slopes  as  nodal  degrees 
of  freedom.  By  selecting  an  appropriate  value  of  the  penalty  parameter, 
solutions  for  thin  as  well  as  thick  plates  are  obtained.  Numerical  results 
are  presented  for  rectangular  plates  with  various  edge  conditions. 
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1.  INTRODUCTION 

The  classical  thin-plate  theory  assumes  that  normals  to  the  midsurface 
before  deformation  remain  straight  and  normal  to  the  midsurface  after  de- 
formation, implying  that  the  transverse  shear  deformations  are  negligible. 

As  a result  the  free  vibration  frequencies  calculated  using  the  thin-plate 
theory  are  higher  than  those  obtained  by  the  Mindlin  plate  theory  [1],  which 
includes  transverse  shear  flexibility  and  rotary  inertia  effects;  the  devia- 
tion increases  with  increasing  mode  number.  The  transverse  shear  effects 
are  even  more  pronounced,  due  to  the  low  transverse  shear  moduli  relative 
to  the  in-plane  Young's  moduli,  in  the  case  of  filamentary  composite  plates. 
A reliable  prediction  of  the  response  characteristics  of  high-modulus  com- 
posite plates  requires  the  use  of  shear  deformable  theory. 

A number  of  shear  deformable  theories  have  been  proposed  to  date.  The 
first  such  theory  for  laminated  isotropic  plates  is  apparently  due  to  Stavsky 
[2].  The  theory  has  been  generalized  to  laminated  anisotropic  plates  by 
Yang,  Norris  and  Stavsky  [3].  A good  review  of  various  other  theories  (e.g., 
[4,5])  can  be  found  in  [6].  It  has  been  shown  (see,  for  example,  [5-8])  that 
the  Yang-Norris-Stavsky  (YNS)  theory  is  adequate  for  predicting  the  flexural 
vibration  response  of  laminated  anisotropic  plates  in  the  first  few  modes. 
Whitney  and  Pagano  [9]  applied  the  YNS  theory  to  the  cylindrical  bending 
of  antisymmetric  cross-ply  and  angle-ply  plate  strips  under  sinusoidal 
lading  and  free  vibration  of  antisymmetric  angle-ply  plate  strips  (see  also 
[10,11].  More  recently,  Bert  and  Chen  [12]  presented,  using  the  YNS  theory, 
a closed-form  solution  for  the  free  vibration  of  simply  supported  rectangular 
plates  of  antisymmetric  angle-ply  laminates. 

While  considerable  effort  has  been  expended  in  the  finite-element  vibra- 
tion analysis  of  isotropic  plates,  only  limited  investigations  of  laminated 
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anisotropic  plates  can  be  found  in  the  literature  [13-16].  Exploiting  the 
symmetries  exhibited  by  anisotropic  plates,  Noor  and  Mathers  [13,14]  studied 
the  effects  of  shear  deformation  and  anisotropy  on  the  accuracy  and  convergence 
of  several  shear-flexible  displacement  finite-element  models  based  on  a form 
of  Reissner's  plate  theory.  The  analysis  was  limited  to  symmetrically  lami- 
nated cross-ply  plates.  In  [15,16]  vibration  of  only  cross-ply  laminated 
plates  was  considered. 

All  of  the  plate  bending  elements  based  on  the  displacement  type  formu- 
lations are  algebraically  complex  (require  (^-continuity;  see  Kawai  and 
Yoshimura  [22])  and  are  computationally  too  expensive  to  use  in  nonlinear 
analyses.  Several  C°  elements  have  been  developed  in  the  last  decade  in 
the  interest  of  computational  efficiency  in  large  nonlinear  problems.  These 
inlcude  the  so-called  mixed  finite  elements  (see,  Hermann  [23],  Nemat-Nasser 
and  Lee  [24],  Kikuchi  and  Ando  [25],  and  Reddy  and  Tsay  [26]),  and  elements 
based  on  the  discrete  "Kirchhoff  hypothesis"  of  Wempner,  Oden  and  Kross 
[27]  and  Fried  [28],  and  on  the  "reduced  integration"  techniques  of  Fried 
[29],  Zienkiewicz  and  Hinton  [30],  and  Hughes,  Cohen,  and  Haroun  [31].  In 
these  works  only  linear  plate  bending  problems  were  considered. 

The  purpose  of  the  present  investigation  is  two-fold:  to  develop  a finite- 
element  model  based  on  Yang-Norris-Stavsky  (YNS)  theory,  and  to  present  a 
penalty-finite  element  that  treats  the  slope-displacement  relations  as  con- 
straints. The  relationship  between  the  penalty  formulation  and  Reissner- 
Mindlin  thick-plate  equations  is  established.  Numerical  results  are  presented 
for  free  vibration  of  antisymmetric,  angle-ply  laminated-composite  plates 
(based  on  YNS  theory),  and  for  nonlinear  bending  (based  on  von  Karman  theory) 
of  thin  plates.  Results  are  compared  with  the  closed-form  and  other  approxi- 
mate solutions. 
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Following  this  introduction,  we  review  the  equations  of  the  Yang-Norris- 
Stavsky  theory,  and  present  an  associated  variational  formulation,  which  is 
required  for  the  finite -element  modeling. 

2.  GOVERNING  EQUATIONS  AND  VARIATIONAL  FORMULATION 
Consider  a plate  of  constant  thickness  h composed  of  an  even  number 
of  thin  anisotropic  layers  oriented  alternately  at  angles  6 and  -e.  The 
origin  of  the  coordinate  system  is  located  at  the  middle  plane  of  the  plate 
with  the  z-axis  being  normal  to  the  mid-plane.  The  material  of  each  layer 
is  assumed  to  possess  a plane  of  elastic  symmetry  parallel  to  the  xy-plane. 
We  shall  denote  the  middle  plane  with  .1. 

The  YNS  theory  is  based  on  the  following  assumed  displacement  field: 


u = uQ(x,y,t)  + z vx(x,y,t) 
v = v0(x,y,t)  + z yy(x,y,t) 


0 


w = w(x,y,E) 

where  u,  v,  and  w are  the  displacement  components  in  the  x,  y and  z-d'i  recti ons , 

respectively,  uQ  and  vQ  are  the  in-plane  (stretching)  displacements  of  the 

middle  plane,  and  q and  v are  the  shear  rotations.  Recalling  the  strain- 
x y 

displacement  equations  of  linear  elasticity,  we  have 
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Owing  to  the  existence  of  a plane  of  elastic  symmetry,  the  constitutive 
relations  for  any  laser  in  the  (\,y'  system  are  given  by 
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ivnere  • are  the  plane-stress  reduced  stiffness  components  of  the  layer  material 
Introducing  the  stress  and  moment  resultants  per  unit  length. 
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With  the  finite-element  discretization  in  mind,  we  give  a variational 
form  of  the  equations  in  (7).  The  kinetic  energy  and  the  strain  energy 
associated  with  the  problem  are, 

T = I j ^ + (]f)  3 + ICtjX*)  + ]*  dA  dt 

U * i f ft/2  j • . dA  dz  , 

2 Jfl  J 'J  U 
* -t/2 

Using  equations  (2)  and  (4),  one  can  express  the  strain  energy  in  terms  of 
the  displacements  and  the  slope  functions  as 
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It  can  be  verified  that  the  first  variation  of  the  Lagrangian  functional,  L=T-V 
(i.e.  Hamilton's  principle)  leads  to  the  equations  of  motion  (7)  expressed 
in  terms  of  the  displacements  and  slope  functions.  Here  V denotes  the  total 
potential  energy  (i.e.  sum  of  the  strain  energy  and  energy  due  to  applied 
loads)  of  the  plate. 


3.  FINITE -ELEMENT  MODELS  OF  THE  YNS  THEORY 


The  first  variation  of  the  Lagrangian  functional,  L = T-V,  gives  the 
variational  form  of  the  equations  of  motion  governing  the  layered  aniso- 
tropic plates.  This  variational  form  is  convenient  for  the  finite-element 
formulation.  Since  we  are  primarily  interested  here  in  the  free  vibration 
analysis,  the  potential  energy  due  to  the  applied  loads  is  zero. 

Let  the  domain  Q.  be  decomposed  into  a set  of  finite  elements.  The 
restriction  of  the  Lagrangian  functional  L to  the  finite  element  a is 
denoted  by  L . We  have 


L 


N 

= Z 
e=l 
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where  N denotes  the  total  number  of  finite  elements  in  the  mesh.  Over  each 
element  flg,  the  field  variables  (uQe  , voe  , we,  ^xe  , i e ) can  be  represented 
by  the  following  approximation 
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where  (x,y)  are  finite-element  interpolation  functions  (also  cal  led  shape 
functions).  For  simplicity,  the  same  interpolation  is  employed  for  all  of 
the  five  fields.  Substitution  of  equation  (12)  into  the  first  variation 
.1 , of  L gives  (omitting  the  algebraic  details)  the  following  set  of  equa- 
tions for  a tvpical  element  : : 
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The  elements  (x,-:  = 1,2...., 5 ; i ,j  = l,2....,r)  of  the  stiffness  matrix 
and  m2,  of  the  mass  matrix  are  defined  by 
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.Note  that  in  the  above  derivation,  the  parameter  A^  is  set  to  cero  for 
s i mp 1 i c i ty . 

For  free  vibration,  the  (discrete)  equations  of  motion  (1 3) , after  assem- 
bly of  the  elements,  take  the  form  of  a standard  eigenvalue  problem: 

([<]  - ."[M])Ul  = i0'  (16) 

where  is  the  frequency  of  the  natural  vibration.  Equation  (16) 
can  be  solved,  after  imposing  the  boundary  conditions  of  the  problem,  by  any 
standard  eigenvalue  program. 

Several  comments  are  in  order  on  the  finite-element  equations  (16).  In 
theory,  each  layer  of  the  plate  can  have  its  own  material  properties  and 
angle  of  orientation.  The  material  properties  of  each  layer  are  generally 
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known  in  its  own  coordinate  system.  Hence,  they  must  be  transformed,  using 

the  angle  of  orientation,  to  the  plate  coordinates  (see  [17])  to  get  Qy^ 

for  each  lamina . To  obtain  A . • , B . ■ , and  D • . one  must  use  equations  (10). 

* J 1 J • J 

For  example. 
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where  Z„,  is  the  distance  from  the  mid-plane,  alonq  the  thickness,  of  the  lower 

ill 

surface  of  the  m th  layer.  Thus,  once  A..,  B..  and  0..  are  known  for  the  lay- 
ered  composite  plate,  it  can  be  viewed  as  a plate  made  of  single  material 
having  behavior  characterized  by  the  constants  Ay,  By  and  Dy . Con- 
sequently, the  finite  element  procedure  becomes  the  one  that  is  used  in 
Mind! in's  plate  theory. 

In  the  present  study  the  eight-node  (r=8)  rectangular  elements  (of  the 
serendipity  family)  are  used.  The  element  stiffness  matrix  is  of  order 
40  by  40. 


4.  FREE  VIBRATION  ANALYSIS  BASED  ON  THE  YNS  THEORY 
Numerical  results  are  presented  for  antisymmetric  angle-ply  rectangular 
plates  with  all  four  edges  simply  supported.  All  of  the  layers  are  assumed 
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In  addition  to  the  isotropic  properties,  the  following  two  sets  of  dimension- 
less material  properties  (typical  of  two  kinds  of  graphi te-epoxy)  are  used: 


1. 

Isotropic: 

Ei/E2  - 

1.0  , 

G]2/E2  = 1/2(1 +v ) , v = 0.3 

2. 

Material  I: 

El/E2  = 

40  , 

G>12^E2  = 0.6  , G-j  g/E2  = G23/E2  = 0.5, 

3. 

Material  II 

: El/Ez- 

25  , 

Gi2/e2  = 0.5  , G-|  3/^2  = 0*2  » v12  = 

All  of  the  computations  were  carried  out  on  an  IBM  370/158  in  double  precision. 
Due  to  storage  limitations,  onlyaquarter  plate  (exploiting  the  biaxial 
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Figure  1 Boundary  conditions  and  the  finite  element  mesh  for 
the  simply  supported  plate. 

symmetry)  is  used  to  compute  the  fundamental  natural  frequencies  and  a half- 
plate  model  is  used  to  obtain  the  symmetric  and  antisymmetric  higher  order 
modes.  The  fini te -element  mesh  and  the  boundary  conditions  for  the  simply 
supported  edges  are  shown  in  Figure  1. 

As  a check  for  the  numerical  accuracy  of  the  finite  element  method, 
natural  frequencies  were  first  obtained  fora  thick  isotropic  plate  (v  = 0.3) 
with  side-to-thickness  ratio  of  10.  The  finite-element  results  are  compared 
with  3-D  linear-elasticity  solution,  Mindlin's  thick-plate  theory,  and  the 
classical  thin-plate  theory  in  Table  1.  The  finite-element  solution  was 
obtained  using  2x2  and  4x2  quadratic  element  meshes  in  the  half  plate. 

Clearly,  the  results  obtained  by  using  the  4x2  mesh  are  closer  to  the  thick- 
plate  theory  indicating  that  further  mesh  refinements  could  yield  better 
accuracy.  However,  with  the  increasing  mode  number, the  finite-element  results 
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Comparison  of  nondimensional  frequencies 


^C(c/Er) 


1/2 


of  a square  simply  supported  isotropic  plate  (v  = 0.3,  a/t  - 10) 
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74.041 
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Table  2 Effects  of  in-plane  displacements,  lamination  anqle,  and  the 

finite-element  mesh  on  the  dimensionless  fundamental  frequency, 

J 9 \ 

\ = ^a"( j/E , t" ) ',  of  a simply  supported  square  plate 
(a/t  = 10,  a/b  =1,4  layers) 


! 
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r 

< 

V 

NDF  = 3 
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NDF- 3 

NDF  = 5 

[12] 
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Ll 
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1x1 

14.504 
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13.492 

19.268 

19.412 
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20. 183 
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9 v9 
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14.193 
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17.63 

19.154 

18.609  18.46 
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T I 

lxl 

10.539 

11.978 

13.164 

13.543 

13.738 

13.993 

L L 

! 

9 v9 
2x2 

10.358 

10.358 

13.013 

12.739 

13.591 

I 1 

13.303) 

[ 
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seem  to  deviate  more  from  the  3-D  elasticity  and  Mindlin's  thick  plate  solu- 
tions . 

Figure  2 shows  a plot  of  the  dimensionless  frequency  versus  the  side-to- 
thickness  ratio  for  three  aspect  ratios  of  layered  composite  plate  consisting 
of  four  Material  I layers  oriented  at  45°/-45°/45c’/-450 . The  results  are 
compared  with  the  closed-form  solution  of  Bert  and  Chen  [12].  For  a given 
aspect  ratio,  the  finite-element  results  seem  to  approach  the  closed- form 
solution  with  a decreasing  side-to- thickness  ratio.  The  figure  also  indicates 
that  the  finite-element  solution  converges  to  the  exact  with  decreasing 
aspect  ratio.  Figure  2 also  contains  a plot  of  the  fundamental  frequency 

■ 

versus  the  side-to-thickness  ratio  of  a square  plate  consisting  of  four- 
layers  of  Material  II.  All  of  the  results  were  obtained  using  a 2x2  mesh 
of  quadratic  elements  in  the  quarter  plate. 

Table  2 shows  the  effect  of  in-plane  displacements,  the  lamination  angle, 
and  the  mesh  on  the  fundamental  natural  frequency  of  square  laminated  plates 
consisting  of  four  layers  (a/t  = 10).  The  three-degree-of- freedom  (NPF  = 3) 
solution  is  obtained  by  suppressing  the  inplane  displacement  degrees  of 
freedom  (u,v).  It  is  clear  that  the  fundamental  frequency  increases  with 
the  lamination  angle  (results  are  symmetric  with  respect  to  45°,  i.e., 
results  for  90°-  > are  identical  to  those  of  o),  and  the  3-degree-of- freedom 
(w,i|i x,>y)  prediction  results  in  higher  frequencies.  Due  to  storage  limita- 
tion further  refinement  of  the  mesh  was  not  possible;  however,  the  finite- 
element  solution  seems  to  converge  toward  the  closed-form  solution. 

The  effects  of  the  lamination  angle  and  the  number  of  layers  on  the 
dimensionless  fundamental  frequency  are  shown  in  Table  3.  The  finite 
element  method  predicts  the  lower  values  for  the  frequencies  compared  to  the  closed 
form  solution  of  Bert  and  Chen  [12]  as  the  number  of  layers  is  increased. 


________ 


J 


0.0  10.0  CO.O  30.0  -10.0  50.0  60.0 


Figure  2 Comparison  of  fundamental  frequency  for  a four-layer  antisymmetric 
angle-ply  ( e =45 3 ) simply  supported  square  plate  by  classical  plate 
theory  (CPT)  and  shear  deformation  theory  of  Ref.  [3]. 

Table  3 Effects  of  lamination  angle  (e)  and  number  of  l^y^rs  on  the 

dimensionless  fundamental  frequency  _a-’(i'/C;t:)  *'  “ of  a simply 
supported  square  plate  ( 3/-6/8/. . ./-c , a/t  = 10) 


; No. 

or 

! Lavers 

i 

Material 

Type  I 

1 

Material  Type  II  i 

30 

0 

45 

0 

30° 

1 

45° 

Present* 

[12] 

Present* 

C i on 

Present* 

2 

15.001 

15. 715 

10.04 

11.351 

n.75;  ; 

4 

17.689 

17.63 

18.609 

13.46 

12.739 

13.303 

, 

6 

IS. 002 

18.23 

13.925 

19.09 

12.893 

13.471  i 

1 3 

18.104 

18.42 

19.028 

19.29 

12.950 

13.524 

i 

10 

18.150 

18.51 

19.074 

19.33 

12.973 

I 

13.549  ! 

I >- 

18.175 

~~ 

19.093 

12.985 

13.562 

! 14 

13.1S9 

— 

19.113 

12.993 

13.569 

! 16 

13.199 

18.60 

19.122 

19.48 

12. 498 

13.575  | 

2x2  quadratic  element  is  used 
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Tables  4 and  5 show  the  effects  of  plate  aspect  ratio  and  length-to- 
tnickness  ratio  on  the  dimensionless  fundamental  frequency  (recall  Figure  2) 
for  Material  I while  Tables  6 and  7 contain  similar  results  for  Material  II. 

The  results  are  obtained  using  2 by  2 mesh  of  quadratic  elements.  From 
Table  4,  it  is  clear  that  the  finite-element  results  are  closer,  but  higher, 
than  the  closed-form  solution  of  Reference  [12], 

Finally,  Table  8 shows  a comparison  of  nondimensional  frequencies,  for  a 
four  layered  laminated  square  plate  ( 45°/-450/450/-45° , a/h  = 10),  obtained 
by  various  investigators.  The  table  includes  the  closed-form  solution  ot_ 

Bert  3nd  Chen  [12],  classical  plate  theory  [20],  and  the  finite-element 
solutions  obtained  using  various  meshes  and  with  and  without  in-plane  dis- 
placement degrees  of  freedom  (DOF).  While  the  4 by  2,  3-OOF  model  gives 
better  results  than  the  2 by  2,  3 DOF  model,  the  2 by  2,  5-DOF  model  aives 
more  accurate  results  when  compared  to  the  closed-form  solution. 

While  the  present  study  concentrated  on  the  antisymmetric,  angle-ply 
laminated  plates  with  simply-supported  edge  conditions,  the  theory  oresented 
(and  the  computer  program  developed)  is  valid  for  general  laminated  plates,  and 
for  various  other  edge  conditions  (see  Reddy  [21]). 


5.  PENALTY  FORMULATION  OF  VON  KARMAN  PLATE  EQUATIONS 
Let  the  middle  plane  f.cR-  of  the  plate  coincide  with  the  (x,y)  plane, 
and  let  jjj  be  its  piecewise  smooth  boundary.  The  strain-energy  functional 
associated  with  the  nonlinear  (von  Karman)  theory  of  thin,  isotropic  plates  is 


U (w)  - * | 

C L J , U- 


(*.Xx>2  + (w’yy);  + 2vw’xxw’yy  + (w,xy):~j  dxdy 


+ U (w,  w,  ) 
n^  ’x  ’y' 


’ 1 ^ \ 
V 1 o ) 


wnere  Un(-)  is  the  contribution  due  to  the  large-deflection  assumption. 


Table  4 Effects  of  plate  aspect  ratio  (a/b)  and  length-  to- thickness  ratio  (a/t)  on  the 

2 2 1 

dimensionless  fundamenta 1 frequency,  >.  = ua  (f/E-t  ) 2 of  a simply-supported 
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Table  5 Effects  of  plate  aspect  ratio  (a/b)  and  length- to- thickness 

ratio  (a/t)  on  the  dimensionless  fundamental  frequency, 

2 2,’ 

\ = u>a  (r/E^t  ) ,J  of  a simply-supported  rectangular  plate  made 
of  Material  I (30<7-30°/307  - 30°). 


a/b 

a/t' 

0.2 

0.4 

0.6 

0.8 

1.0 

1.2 

1.4 

1.6 

1.8 

“I 

2.0 

10 

11.106 

12.172 

13.734 

15.615 

17.689 

19.880 

26.275 

24.460 

26.836 

29.2S8  ! 

I 

20 

12.595 

13.950 

15.969 

18.457 

21.281 

24.362 

27.645 

31 . OSS 

34.6S0 

38 .428  - 

30 

12.972 

14.401 

16.543 

19.204 

22.253 

25.619 

29.244 

33.086 

37.132 

41 .390 

40 

13. 136 

14.591 

16.779 

19.507 

22.646 

26.123 

29.896 

33.909 

38.154 

42.641 

50 

13.233 

14.700 

16.908 

19.667 

22.850 

26 . 3S9 

30.229 

34.330 

38.677 

43.284  ' 

60 

13.301 

14.773 

16.992 

19.767 

22.974 

26.545 

30.426 

34.577 

3S.985 

4 j . fco  j 

Table  7 Effects  of  plate  aspect  ratio  (a/b)  and  length- to- thickness 

ratio  (a/t)  on  the  dimensionless  fundamental  frequency, 

2 2 1 

\ = oia  (a/E9t  ) *,  of  a simply-supported  rectangular  plate  made 
of  Material  II  (30o/-307  30o/-30°). 


L^a/b 

a/t""- 

0.2 

0.4 

0.6  0.8 

: 1.0 

1.2 

1.4 
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1 .S 
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CO 
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! 60 
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10.555 

11.731 

13.505 

15.728 
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Table  C Effects  of  plate  aspect  ratio  (a/b)  and  length- to- thickness  ratio  'a/t)  on 

2 2 2 

the  dimensionless  fundamental  frequency,  • = ..a  (. /E^t  )\  of  a simply-supported 
rectangular  plate  made  of  Material  II  (45''/-45°/45'7-45c' ) 
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the  thickness  or'  me  elate.  Introducing  the  slopes  •>  and  , , 


" - > • 0 
■ \ \ 


- A 

\\  V 
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U (d  ,8  ) = 5 f 

8 ' x y 2 j 


a - + p-  + 2v  0 0 + t'—  (o  + 0 )—  I dxdv 

X.X  y,y  x,xLy,y  2 v x,y  y,x;  ux  - 


+ U (0  6 ) 

n ' x , y ' 


(22 ) 


Variation  of  the  functional  in  (18)  leads  to  the  equilibrium  equation  in 
terms  of  the  deflection  w,  while  that  of  (22)  leads  to  equilibrium 
equations  in  terms  of  the  slopes.  The  strain  energy  associated  with  the 
thick  plate  (Reissner's)  theory  is  given  by 


x * vy(y  * ‘-V>'x,xvy ,y  ^ 2 ^^x^v  ^ , x ^ 


dxdy 


+ ktG  1 


7 i [(w.x  - *x)-'  f (w,  - ,v);]  dxdy  + Un(;x,v  ) 


(23) 


where  ,x  and  vy  are  the  slope  functions  and  k is  the  shear  correction  factor. 

Note  the  similarity  of  the  expressions  in  (22)  and  (23). 

In  the  penalty  method  the  problem  of  finding  the  critical  points  (e  ,0  ) 

x y 

of  the  functional  U subject  to  the  constraint  condi tions (21 ) is  formu- 

U 

lated  as  one  of  seeking  the  critical  points  (ox,o  ) of  the  modified  function- 
al. 


up(6x.0y.w)  = U0(0x,0y) 


Jo 


+ (w,, 


,)•']  dxdy 


24) 


without  constraints.  Here  t-  is  the  penalty  parameter.  The  penalty  para- 
meter c plays  a crucial  role  in  the  accuracy  of  the  method.  In  theory,  as 
r goes  to  infinity,  the  solution  (eyU),  0..(e))  to  the  penalty  problem 
approaches  the  true  solution  (ox,  o ) (see  Reddy  [33]).  The  shear  force 
resultants  can  be  identified  as 


K - 


ex(c)) 


(25) 


, Qj  - c(M^-  ey(0) 


Also  note  that  the  penalty  functional  Uq  is  of  the  same  form  as  the  function- 
al Us  associated  with  shear  flexible  theory  with  the  following  correspondence: 


'i'x  — 


y 


- ktG 


(26) 


Thus,  for  a very  large  value  of  t we  recover  the  strain  energy  of 
the  thin-plate  theory  from  equation  £4),  and  for  t = ktG,  represents 
the  strain  energy  associated  with  Reissner’s  thick- plate  theory. 


6.  PENALTY  FINITE-ELEMENT  MODEL 

The  total  potential  energy  associated  with  the  penalty  formulation  of 
the  equations  governing  the  large  deflection  of  elastic  (isotropic)  plates 
is  given  by 


Tp(w,9x,9y)  = up<ex>Vw>  * i . wP  dxdy  ' j iwds  " i Mn  Tn  ds 

0 3i'-m 


where  3fl  + = 3L,  P is  the  distributed  load,  and  n is  the  unit  outward 

normal . 

We  seek  the  finite-element  solutions  of  the  form 


w = iw.n.  , e = se’N-  , e>(  = eq’n. 
it  x x i ’ y y i 


(28) 


where  are  the  (finite  element)  shape  functions.  Substituting  (2S)  into 
(27),  we  obtain  the  following  equations  for  a typical  element: 


[Kx]  + e[s'] 


CKXy] 


l u 


' V 


Q) 


(29) 


where 
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[<x]  = 0([SX]  ♦ (l-v)[Sy])  + [Gx]  , Sf]  = 


N,  _N  • dvdn 
‘ > Jo 


[Ky]  - D([Sy]  ♦ (l-v)[Sxj)  * (Gy]  , [ Kxy]  = D[Sxy]  + [Gxy] 


Gx  . . f 

J 


, vvj dxdy  * g:!j  a viNj dxdy  ’ G?j s ; 

*•  * ft 


Vi"j  d*dy 


[K]  ■ [s’1]  ♦ [Sy]  . B*  = ' «x»,ds  , B?  . ’ 


..  Vids 


Computational  Aspects.  Several  comments  are  in  order  on  the  computa- 
tional aspects  of  the  penalty  finite-element  model  developed  herein.  First 
let  us  examine  the  consequence  of  the  limit  £ 1 going  to  zero.  The  penalty 
finite-element  equation  (12)  can  be  expressed  in  the  form 

l [ 3 + £ [ Ko  ] ) '.  d : = { F • 

For  very  large  values  of  £ , the  second  term  (i.e.  penalty  term)  in  the 
parentheses  dominates  and  we  have, 

7 --F;  = ~ ([K,  ] + e[K-];d;  : [K7](d;  , and  lim  [K,]{d'>  = \0; 


Thus,  for  a very  large  value  of  e,  the  bending-energy  contribution^  the 
stiffness  matrix  is  lost  (and  the  energy  due  to  transverse  shear  dominates), 
and  equation  (32)  yields,  in  an  attempt  to  satisfy  the  constraint  conditions 
exactly,  the  trivial  solution  (d:  = (0)  (if  [1U]  is  nonsingular).  If  the 
equations  in  (15)  are  linearly  dependent  (i.e.  [\?]  is  singular),  then  the 
number  of  independent  (constraint)  equations  are  less  than  the  number  of 
unknowns  and  the  remaining  equations  are  provided  by  [K-j].  Thus,  it  is 
necessary  to  have  nonsingular  [Ki]  and  singular  [K~]  in  order  to  obtain  a 
nontrivial  solution.  The  same  conclusion  can  be  reached  from  a mathematical 
point  of  view.  The  existence  and  uniaueness  of  solutions  to  any  finite 
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element  model  requires  the  associated  bilinear  form  to  satisfy  certain  con- 
tinuity and  coercivity  conditions  (see  Oden  and  Reddy  [34]).  For  the  pen- 
alty finite-element  model,  these  conditions  impose  restrictions  on  the  choice 
of  interpolation  used  for  w,  and  v , and  e to  make  the  parameter  appearing 
in  the  coercivity  condition  be  independent  of  the  mesn  sice,  n.  The  latter 
can  also  be  achieved  by  employing  reduced  integration  on  the  penalty  (or 
shear)  terms,  which  makes  the  matrix  [K?]  singular.  Reduced  integration  is 
employed  in  the  present  study  to  evaluate  the  elements  of  [K-].  For  example, 
if  a bilinear  quadrilateral  element  is  used,  a lxl  Gauss  rule  must  be  used 
for  [K,]  in  place  of  the  standard  2x2  Gauss  rule.  As  for  the  value  of  . , 
it  snculd  be  equal  to  kGt  for  thick  plates  (in  which  i - ,x  t sw/ix,  etc.), 
and  an  order  or  two  greater  for  the  thin-plate  analysis  (to  satisfy  the 
constraints  Jx  = ?w/<x,  and  = i-w/.-y). 

7.  NONLINEAR  BENDING  AND  FREE  VIBRATION  RESULTS 
The  penalty  finite  element  described  above  is  used  in  the  (static) 
oencing  analysis  of  square  plates  under  uniformly  distributed  load  for 
various  edge  conditions.  In  these  studies  a value  of  0.3  is  used  for  tne 
Poisson's  ratio,  and  c is  taken  to  be  Gt/n,  where  h is  the  element  (mesh) 
sice.  Due  to  the  biaxial  symmetry  of  the  proolems  considered  here,  only 
quarter  plates  are  anaiyced  in  the  interest  of  computational  efficiency  and 
convenience.  The  following  iterative  scheme  is  used  to  obtain  the  con- 
vergence of  the  transverse  deflections.  Beginning  with  the  cero  solution, 
tne  linear  solution  is  obtained  after  the  first  iteration.  This  solution 
is  used  to  compute  the  nonlinear  (geometric)  stiffness  coefficients  for  the 
next  iteration.  Only  a weighted  average  of  the  following  form  is  used  to 
calculate  the  geometric  stiffness  coefficients: 


(33 
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where  id  is  the  solution  vector,  n is  the  iteration  number,  and  a is  a 
parameter,  0 < * < 1.  Preliminary  investigation  showed  that  a value  of 
= 2/3  gives,  in  the  present  formulation,  the  fastest  convergence.  The 
convergence  criterion  used  is  of  the  form, 

[ E (aJ+1  - a")-]:/:  < e = 10'4  (34) 

All  of  the  computations  are  carried  on  an  IBM  370/158  model  in  double  pre- 
cision. Due  to  the  space  limitation  only  limited  results  are  included  here. 

Figures  3 and  4 show  the  plots  of  the  (normalized)  center  deflection 
for  simply-supported  and  clamped  plates,  respectively,  under  uniformly 
distributed  loading,  and  corresponding  stresses  are  shown  in  Figures  5 and 
6.  The  present  solutions  are  in  good  agreement  with  those  of  Way  [35], 

Wang  [36],  Levy  [37],  and  Yamaki  [38].  The  present  solutions  (for  various 
meshes  and  elements)  are  bounded  by  Yamaki 's  solution  from  above,  and 
Wang's  solution  from  below  for  simply-supported  plate. 

Table  9 shows  a comparison  of  the  penalty  (present)  solution  with 
the  mixed  solution  of  Reddy  and  Stricklin  [39]  for  the  clamped-hinged  square 
plate  under  uniformly  distributed  load.  The  mixed  finite  element  in  [39] 
was  based  on  the  von  Karman  equations.  The  bending  moments  (M  = My) 
are  compared,  for  the  same  problem,  in  Table  10.  The  values  of  the  bending 
moment  are  obtained  at  the  nodes  in  the  case  of  the  mixed  method,  whereas 
they  are  computed  at  the  Gaussian  points  in  the  case  of  the  penalty  method. 

The  locations  of  the  Gaussian  points  are  given  in  parentheses. 

Due  to  the  simplicity  of  the  present  element,  it  is  very  economical 
to  use  in  nonlinear  problems  which  require  large  numbers  of  elements.  The 
complicated  conventional  bending  elements  (i.e.  elements  based  on  displacement 
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Figure  4 Load-deflection  curves  for  clamped  square  plate  under 

uniform  loading 
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Figure  5 Stress-load  curves  for  simply  supported  square  pla 

under  uni  to  rm  loading 
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Figure  6 Stress-load  curves  for  clamped  square  plate  under 
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type  formulations)  are 
tions  are  compared  but 


only  marginally  more 
computati onal ly  very 


accurate  when 
expensi ve. 


the  1 i near  sol u- 


Table  9.  Comparison  of  Center  Deflection  (w/t)  Obtainej  by  the 
Mixed  and  Penalty  Methods  for  Cl amped-Hi nged  Square 
Plate  Under  Uniformly  Distributed  Load  ( =0.3,  a/t  = 10) 


[ 

LOAD 

1 

MIXED  FEM  [39] 

! 

PENALTY  FEM 

LINEAR 

QUADRATIC 

LINEAR 

QUADRATIC 

Pa-  Et" 

Jx'J 

♦ X 4 

.x^ 

4x4 

oxo 

4x4 

25.0 

i 0 

.494 

0 

451 

0 

436 

0 

.437  , 

0.440 

0.441 

50 . 0 

1 0 

.804 

0 

,724 

0 

, 697 

0 

.699  ■ 

0.707 

0.707 

75.0 

1 

.021 

0 

,916 

0 

,881 

0 

,882 

0.895 

0.393 

100.0  j 

1 

.190  . 

1 

, 066 

1 

,025  ' 

1 

.026  1 

1.043 

1 .040 

125.0 

1 

.330 

1 

1 92 

1 

, 146 

1 

.145  j 
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1.163 

150.0  | 

i 

. 450 

1 

,301 

1 

,252 

1 

,249  j 
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1.270 

175.0 

1 i 

.557 

1 

,398 

1 

.345 

1 1 

.341  j 

1 .371 

1 .365 

200.0  j 

; i 

.652 

1 

,486  1 

1 

.431 

1 

.425  ; 

1 .45S 

1 .452 

225.0 

i i 

.741  | 

1 

. 566 

1 

.509 

1 

, 501  | 

1.539 

1.532 

250.0  | 

l 1 

! i 

i 

.321 

1 

,641 

1 

,581 

1 
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1 .613 

— 

Table  10.  Comparison  of  Maximum  Bending  Moment  (-M  x 10*:)  Obtained 
by  the  Mixed  and  Penalty  Methods  for  Cl amped-Hi need 
Under  Uniformly  Distributed  Load v 1 0.3,  a t s 1 1) 


LOAD 

1 PaVEt4 

MIXED  F 

EM  [39] 

PENALTY 

FEM 

LINEAR 

QUADRATIC 

LINEAR 

6x6 

(.4583,-4533 

; QUADRATIC 

j 4x4 

{(.4375, .4375) 

2x2 

1 

4x4 

2x2 

4x4 

I 25.0 

1 .450 

1.535 

1.617 

1.572 

1 .094 
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2.400 
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2.829 

2 “*?q 

1.327 

1 .420 
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3.553  | 
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2.381 

1 .810 

100.0 
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. 2. IIS 
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5.043 
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5.300 

3.234 

2.375 

150.0 

4.453 

5 . 683 

6.256 

6.023 

3.585 

2.596 
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4.306 

6.275 

6.942 
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3.903 

? 792 
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5.117 

6.327  < 
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4.194 

2.967 
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5.407 

7 . 346 

3.199 

7.986 

4.464 

3.126 
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5.671 

7 . S36 

3 . 777 

3.592 
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1 

3.  APPLICATION  TO  PLATES  MADE  OF  B I MODULUS  MATERIALS 


The  finite  element  model  developed  herein  con  also  be  applied  to  the 
analysis  of  plates  constructed  of  bimodulus  composite  material.  Due  to 
different  elastic  properties  in  tension  (T)  and  compression  (C),  the 
material  coefficients  A.j,  B.j,  and  D.j  are  now  defined  by  (see  [40]), 


(a.  . ,b . . ,d . 

ij  iJ  ij 


•t  12 


— t / 2 


Qij  ( 1 , z , z ■ ) dz  , (i  ,j  = l ,2,6) 


r 'n  .-“t+l 

£ ! Q.  ,dz  + 1 Q.  d; 

s-i  K,  u,t  ci u ’j2i 


(35) 


Here  Q^.  denote  the  stiffness  coefficients  in  the  plate  coordinates  of 
the  v-th  layer  in  tension  (k=l)  or  compression  (k=2),  and  z is  the  distance 
from  the  midsurface  to  the  neutral  plane  (which  is  unknown  a priori). 

Figure  7 shows  the  influence  of  the  aspect  ratio  (b/a)  and  side-to- 
thickness  ratio  (a/t)  on  the  location  of  neutral  surfaces  for  a single-layer, 
isotropic,  bimodulus,  simply  supported  rectangular  plate  subjected  to 
sinusoidal  loading. 


p = Pq  sin(:x/a)  sin(iy/b) 


The  following  elastic  properties  are  used: 

En  = 3.584  GPa  , E^  = 1.792  GPa  , Fr:  = E,i  , E*;  = Eln 


(36) 


Gi:  = d\2  = 1 -27  GPa  . v;:  = = 0.4 


vi2  “ vai  = u.z 

Note  that  for  b/a  = 1,  the  neutral  surfaces  associated  with  x-  and  y-directions 


coincide  (i .e. , 
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Neutral -surface  location  vs.  plate  aspect  ratio 
and  side-to-thickness  ratio  for  single-layered 
rectangular  plate  under  sinusoidal  loading. 


Two-layer  cross-plv 


Neutral -surface  location  vs.  plate  aspect  rati 
for  two-layer,  cross-ply  vO''  90°)  square  plate 
under  sinusoidal  loading. 
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Similar  results  are  presented  in  Figures  8 and  9 for  a two-layer, 
cross-ply  (0"/90L')  rectangular  plate  under  sinusoidal  loading.  The  bi- 
modulus properties  used  are  the  same  as  those  listed  in  Table  3 of  [40]. 

Note  from  Fig.  8 that  the  neutral -surface  location,  znx  is  virtually  un- 
changed for  aspect  ratio  greater  than  1,  while  the  neutral -surface  location, 
zny,  increases  in  proportion  to  aspect  ratio.  It  should  also  be  noted 
that  the  neutral  surfaces  do  not  coincide  in  the  cross-ply  case  for  b/a  = 1. 

Figure  10  shows  the  influence  of  the  aspect  ratio  and  side-to- 
thickness  ratio  on  the  transverse  deflection  for  single-layer  and  two- 
layer  cross-ply  problems  discussed  above.  The  effect  of  thickness  on  the 
deflection  is  more  pronounced  than  the  effect  of  the  aspect  ratio. 

9.  SUMMARY  AND  CONCLUSIONS 

Using  the  penalty  function  concept  of  Courant,  a simple  finite  element 
for  the  Yang-Norri s-Stavsky  (YNS)  theory  of  composite  plates  is  developed. 

The  idea  is  also  extended  to  the  nonlinear  theory  of  thin  elastic  plates 
(due  to  von  Karman).  In  these  formulations  the  slope-deflection  relations 
are  treated  as  constraints.  Rectangular  isoparametric  elements  of  the 
"serendipity"  family  are  developed.  Numerical  results  of  natural  frequencies 
are  presented  for  rectangular  plates  of  antisymmetric  angle-ply  laminates. 

In  the  nonlinear  analysis,  only  thin,  isotropic,  rectangular  plates  are  con- 
sidered in  order  to  compare  the  accuracy  of  the  present  element  with  other 
approximate  solutions.  Finally,  the  penalty  finite  element  based  on  the  YNS 
theory  is  employed  to  study  the  bending  of  laminated,  anisotropic  bimodulus- 
material  plates. 

In  the  study  of  free  vibrations,  the  present  element  gives  very  accurate 
results  for  natural  frequencies.  The  results  are  verified  against  closed- 
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form  solutions.  In  the  case  of  nonlinear  bending  of  thin  plates,  the  results 
are  compared  with  those  obtained  by  a mixed  finite  element  method,  and  other 
approximate  solutions.  The  present  solutions  are  in  close  agreement  with 
those  obtained  by  the  mixed  method,  and  are  bounded  (above  and  below)  bv 
other  approximate  solutions.  For  the  bimodulus  plates,  it  is  shown  that 
the  neutral -surface  location,  even  for  single-layer  plates,  may  vary  consid- 
erably from  the  geometric  midplane,  depending  upon  the  degree  of  bimodularity. 
Also,  the  plate  deflection  is  significantly  affected  by  the  bimodulus  action. 

The  penalty  element  developed  herein  is  computationally  simple  and 
accurate,  and  saves  large  amounts  of  computing  time  in  nonlinear  analyses 
which  require  iteration.  It  is  straight  forward  to  extend  the  present  element 
to  a combination  of  von  karma n and  YNS  theory.  Obviously,  the  element  can 
also  be  used  in  the  analysis  of  shallow  shells.  It  is  the  purpose  of  this 
investigation  to  employ  the  present  element  in  the  analyses  of  plates  and 
shells  constructed  of  bimodulus-material  layers. 
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v material  properties  that  are  typical  of  advanced  fiber-reinforced  composites 
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